#include <QtGui>

#include <iostream>
#include <math.h>

#include "RenderingThread.h"

// use only iteration to derive color
double fromIteration(int iter, double /* zx */, double /* zy */) {
	return iter;
}

// use real part of an orbit
double fromRealPart(int iter, double zx, double /* zy */) {
	return iter+10.0*zx;
}

// use imaginary part of an orbit
double fromImaginaryPart(int iter, double /* zx */, double zy) {
	return iter+10.0*zy;
}

// use sum of real and imaginary part of an orbit
double fromSum(int iter, double zx, double zy) {
	return iter+10.0*(zx+zy);
}

// use angle of an orbit; the end angle between orbit and real/imaginary axis of Gaussian plane
double fromAngle(int iter, double zx, double zy) {
	return iter+5.0*atan2(zx, zy);
}

// more continuous color scheme
// From: Coloring Dynamical Systems in the Complex Plane
// http://math.unipa.it/~grim/Jbarrallo.PDF
double fromNICA(int iter, double zx, double zy) {
	return iter+1-(log(log(sqrt(zx*zx+zy*zy))))/log(2);
}

// it's called by the main UI thread
RenderingThread::RenderingThread(QObject *parent)
: QThread(parent)
{
    mRestart = false;
    mAbort = false;
}

// it's called by the main UI thread
RenderingThread::~RenderingThread()
{
    mutex.lock();
    setAbort(true);
	// wake up the thread if it's sleeping
    condition.wakeOne();
    mutex.unlock();
	
	// we wait until the run method has finished
    wait();
}

// this is the only public method
// it's called by the main UI thread every time a new image is required
void RenderingThread::render(const TRenderParams &renderParams)
{
    // simplyfies locking of a mutex
	// locks when created, unlocks when destroyed (here once the function is left)
	QMutexLocker locker(&mutex);
	
	// copy rendering parameters to this current instance
    this->renderParams = renderParams;
	
    if (!isRunning()) {
		// first time, the thread has to be started
        start(LowPriority);
    } else {
		// thread is either computing new picture right now or it's sleeping, in either cases this will restart it
        setRestart(true);
        condition.wakeOne();
    }
}

// running in it's own thread
// this function never exits until the stop of the thread
void RenderingThread::run()
{
    // the same as for (;;) or Q_FOREVER
	forever {
		// protect accesses to the member variables
		// duplicate all necessary variables in order to minimize the ammount of code needed for locking
        mutex.lock();
		const TRenderParams &renderParams = this->renderParams;
        mutex.unlock();
		
		// let's generate the requested image
		QImage image = renderMandelbrotSet(renderParams);
		if (renderParams.antialiasing)	 {
			// apply filter
			image = medianFilter(&image);
		}
			
		if (aborted()) {
			// delete the renderer
//			delete renderParams.renderer;
			return ;
		}
        if (!restarted()) {
			emit renderedImage(image, renderParams);
		} else {
			// delete the renderer
//			delete renderParams.renderer;
		}
		
        mutex.lock();
        if (!restarted()) {
			// put thread into sleep
			// it will be woken up on a request by the render function
            condition.wait(&mutex);
		}
        setRestart(false);
        mutex.unlock();
    }
}

#define ARGB(qrgb) ((qAlpha((qrgb)) & 0xFF) << 24) | ((qBlue((qrgb)) & 0xFF) << 16) | ((qGreen((qrgb)) & 0xFF) << 8) | ((qBlue((qrgb)) & 0xFF) << 0)
#define QRGB(argb) qRgba(((argb) >> 16) & 0xFF, ((argb) >> 8) & 0xFF, (argb) & 0xFF, ((argb) >> 24) & 0xFF)

// called from it's own thread
// it's actually bit easier to compute each part of imaginary numbers separate than to use std::complex
inline QImage RenderingThread::renderMandelbrotSet(const TRenderParams &renderParams)
{
	QSize size = renderParams.size;
	double scale = renderParams.scale;
		
	// create a new image on the stack
	QImage image(size, QImage::Format_RGB32);
	
	// denotation of which part of the complex plane will we be calculating the set
	double xmin = (renderParams.centerx - size.rwidth() / 2.0) / scale;
    double ymin = (renderParams.centery - size.rheight() / 2.0) / scale;	
    double xmax = (renderParams.centerx + size.rwidth()  / 2.0) / scale;
    double ymax = (renderParams.centery + size.rheight() / 2.0) / scale;
	
//	std::cout << renderParams.centerx << "," << renderParams.centery << "," << scale << ": " << " " << xmin << "," << ymin << "," << xmax << "," << ymax << std::endl;
	
	double zx, zy, zx2, zy2;                    
    double cx, cy;                              
    double cx0, cy0;
                                                
    cy0 = ymin;

	// for all rows in the image
    for (int y = 0; y < size.rheight(); y++) {           
        cx0 = xmin;
		// this get's the pointer for the exact line in the image which can be directly accessed to set a color. It's like big array of colors
		uint *scanLine = reinterpret_cast<uint *>(image.scanLine(y));
		// for each pixel in a row
        for (int x = 0; x < size.rwidth(); x++) {
			// now we start computing the point whether it's a part of the set
            cx = cx0;                             
            cy = cy0;
            zx = zy = 0.0;
			int i;
			
            for (i=0; i < renderParams.maxiter; i++) {
				// periodically check whether the request hasn't been cancelled
				if (restarted() || aborted()) {
					return image;
				}

				// z^2
                zx2 = zx * zx;
                zy2 = zy * zy;

				// check divergation
                if (zx2 + zy2 > 4.0) 
					break;
				// z(n+1)
                zy = 2.0 * zx * zy + cy;
                zx = zx2 - zy2 + cx;				
            }
			
			// draw computed pixel pixel
			QRgb rgb;
			if (i != renderParams.maxiter) {
				rgb = renderParams.renderer->render(i, zx, zy);
            } else {
				rgb = Qt::black;
			}

			*scanLine++ = rgb;
			
			// new point int a line
            cx0 += (xmax - xmin) / size.rwidth();
        }	
		// new line
		cy0+=(ymax-ymin) / size.rheight();
	}
	
	return image;
}

// a very simple implementation of median filter
// called from it's own thread
inline QImage RenderingThread::medianFilter(const QImage *original) 
{	
	QImage image(original->size(), QImage::Format_RGB32);	
	uint *bits = reinterpret_cast<uint *>(image.bits());
	int height = image.size().rheight();
	int width = image.size().rwidth();

	// compute median of each canal separately
	TFilterMatrix r;
	TFilterMatrix g;
	TFilterMatrix b;

	// creates a matrix of each pixel's 3x3 surroundings
	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			int k = 0;
			for (int dy = -1; dy <= 1; dy++) {
				int iy = y+dy;
				if (0 <= iy && iy < height) {
					for (int dx = -1; dx <= 1; dx++) {
						int ix = x+dx;
						if (0 <= ix && ix < width) {
							QRgb rgb = original->pixel(ix, iy);
							r[k] = qRed(rgb);
							g[k] = qGreen(rgb);
							b[k] = qBlue(rgb);
							k++;
						}
					}
				}
			}
			while (k < cFilterMatrixSize) {
				r[k] = g[k] = b[k] = 0;
				k++;
			}
			*bits++ = qRgb(median(r), median(g), median(b));
			// check whether it hasn't been canceled
			if (restarted() || aborted()) {
				return image;
			}

		}
	}
	return image;
}

// The optimized function to find median in a 3x3 integer matrix
// copied from: XILINX XCELL magazine, vol. 23 by John L. Smith
#define PIX_SORT(a,b) { if ((a) > (b)) PIX_SWAP((a), (b)); }
#define PIX_SWAP(a,b) { (a) = (a) + (b); (b) = (a) - (b); (a) = (a) - (b); }
inline int median(TFilterMatrix p) {
	PIX_SORT(p[1], p[2]); PIX_SORT(p[4], p[5]); PIX_SORT(p[7], p[8]);
	PIX_SORT(p[0], p[1]); PIX_SORT(p[3], p[4]); PIX_SORT(p[6], p[7]);
	PIX_SORT(p[1], p[2]); PIX_SORT(p[4], p[5]); PIX_SORT(p[7], p[8]);
	PIX_SORT(p[0], p[3]); PIX_SORT(p[5], p[8]); PIX_SORT(p[4], p[7]);
	PIX_SORT(p[3], p[6]); PIX_SORT(p[1], p[4]); PIX_SORT(p[2], p[5]);
	PIX_SORT(p[4], p[7]); PIX_SORT(p[4], p[2]); PIX_SORT(p[6], p[4]);
	PIX_SORT(p[4], p[2]); 

	return p[4];
}

// --PROPERTIES

inline bool RenderingThread::restarted() const
{
	return mRestart;
}

inline bool RenderingThread::aborted() const
{
	return mAbort;
}
	
inline void RenderingThread::setRestart(bool restart)
{
	mRestart = restart;
}

inline void RenderingThread::setAbort(bool abort)
{
	mAbort = abort;
}